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Method and Apparatus for 
Automated Simulation and Design of Corneal Refractive Procedures 

Field of the Invention 

The present invention relates to systems and techniques for mathematically 
5 modeling a human eye using calculated strain values for a human eye and using a 
mathematical model to simulate strain deformation of the eye by hypothetical incisions, 
excisions, ablations, or prosthetic insertions to arrive at an optimum surgical design by 
identifying the number, shape, location, length, and depth of the incisions, excisions, 
ablations, or of corneal prosthetic insertions required to obtain a uniform or near 
1 0 homogeneous strain pattern on the cornea. 

Background 

The present invention relates to systems and techniques for mathematically 
modeling a human eye using calculated strain values obtained from data measured from a 
human eye. The mathematical model of the present invention simulates the change in 

15 strain conditions of the comea effected by a set of hypothetical incisions, excisions, 
ablations, or corneal prosthetic insertions. A near uniform strain pattern on the comea is 
a critical end-point in the calculation used to arrive at an optimum surgical design for the 
number, shape, location, length, and depth of the incisions, excisions, ablations, or 
corneal prosthetic inserts used in a proposed operation. It should be understood that 

20 hereinafter, including in the claims, the term "incision," which usually refers to a cut 
made by a scalpel, and the term "excision," which usually refers to a cut made by a laser 
beam, are considered to be interchangeable and to have the same meaning. 
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Modem corneal refractive surgery originated with the work of Dr. Svyatoslav 
Fyodorov of Moscow and Dr. Jose Barraquer of Bogota, Columbia. Subsequently, 
various surgical techniques have been developed to alter the curvature of the cornea to 
correct refractive errors. The various techniques include incisional keratotomy using 
5 diamond blades, excisional keratotomy using laser beams to photo-disrupt molecules and 
ablate tissue in a linear pattern, ablative keratectomy or photo-refractive keratectomy 
using laser beams to remove larger areas of corneal tissue, mechanical removal and 
reshaping of corneal tissue (keratomileusis), and implantation of human or synthetic 
materials (corneal prosthetics) into the corneal stroma. All of the known procedures alter 

10 central corneal curvature by changing the structure of the cornea. Additionally, because 
the central corneal curvature is changed, any strain relationships within the comea are 
also changed by these procedures. All such refractive procedures are characterized by 
difficulty in predicting both the immediate and long term results, because of errors in 
calculations of pre-surgical measurements, failure to precisely implement the planned 

15 surgical techniques, and biological variances which affect immediate and long term 
results. 

The comea traditionally has been treated as a spherocylindrical lens, assuming 
that the radius of each individual meridian from the corneal apex to the corneal periphery 
is uniform. Prior methodologies tend to use an approximation to the topographic 
20 information of the comea to determine the refractive power of the eye. In one known 
procedure, circular mires (reflected light images from the comea conventionally used to 
mathematically calculate comeal curvature) are reflected from the corneal surface, and 
the difference between a given point on the mire and an adjacent mire is measured. A 
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semi-quantitative estimate of tlie surface curvature is obtained by comparing this 
measurement with the values obtained using spheres of various radii. Prior mathematical 
models use a variety of approximations such as a simplified form of the corneal surface 
(e.g., spherical) or assume a symmetrical cornea (leading to a quarter model or an 
5 axisymmetric model) or use simplified material properties (e.g., isotropic), or assume 
small deformations or displacements, or do not consider clinically obtained data in the 
construction of the mathematical model. These models, by implicitly assuming uniform 
strain relationships in the cornea, do not accurately model any real strain relationships felt 
by the cornea. 

10 One prior art is the article "On the Computer- Aided and Optimal Design of 

Keratorefractive Surgery," by Steven A. Velinsky and Michael R. Bryant, published in 
Volume 8, page 173 of "Refractive and Corneal Surgery," March/ April 1992. This article 
describes a computer-aided surgical design methodology, proposing that it could be an 
effective surgical design aid for the refractive surgeon, wherein the surgeon could choose 

15 constraints on surgical parameters such as minimum optical zone size, maximum depth of 
cut, etc., measure the patient's corneal topography, refractive error and possible other 
ocular parameters, and then review the computed results. The article refers to several 
mathematical models described in the literature, and how such mathematical models 
might be helpful. However, the article fails to disclose any particular adequate 

20 mathematical model of the corneal strain relationships or any specific recommendation of 
surgical design that has been validated with clinical data. 

Prior keratotomy procedures often are based on the experiential use of 
nomograms indicating appropriate surgical designs for a particular patient based on age, 
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sex, refractive error, and intraocular pressure. These procedures do not account for the 
actual strain relationships in the cornea and frequently result in large amounts of under- 
correction or over-correction. 

Finite Element Analysis (FEA) is a known mathematically based numerical tool 
5 that has been used to solve a variety of problems that are described by partial differential 
or integral equations. This technique has been used primarily in the area of solid 
mechanics, fluid mechanics, heat transfer, electromagnetics, acoustics, and biomechanics, 
including designing remedial techniques being developed for the human eye, to model 
internal structure and stresses in relation to various configurations of intraocular devices 

10 and corneal implants, as described in "Intraocular Lens Design With MSC/pal," by A. D. 
Franzone and V. M. Ghazarian in 1985 at the MSC/NASTRAN User's Conference in 
Pasadena, Calif., and in "Corneal Curvature Change Due to Structural Altemation by 
Radial Keratotomy," by Huang Bisamsin, Schachar, and Black in Volume 110, pages 
249-253, 1988 in the ASME Journal of Biomedical Engineering. Also see "Reduction of 

15 Corneal Astigmatism at Cataract Surgery," by Hall, Campion, Sorenson, and Monthofer, 
Volume 17, pages 407-414, July 1991 in the Journal of Cataract Refractive Surgery. 

There still is an current and continuing need for an improved system for 
accurately predicting outcomes of hypothetical surgical procedures on the cornea to aid in 
the design of minimally invasive corneal surgery. There is a still unmet need for a totally 

20 automated way of determining an optimal design of a surgical plan for incisional, 
excisional, ablative, or insertive keratotomy surgery to meet predetermined visual 
objectives v/ith minimum invasiveness and minimum optical distortion. Further, it would 
be desirable to provide a technique for designing a multi-focal cornea that is similar to a 
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gradient bifocal for patients that have presbyopia. It would be desirable to have an 
accurate mathematical model of the cornea for use in developing new surgical procedures 
without experimenting on live corneas. 

Summary of the Invention 

5 Accordingly, it is an object of the invention to provide a minimally invasive 

surgical procedure for corneal surgery for a human eye to achieve predetermined 
modified characteristics of that eye. 

It is another object of the invention to provide a system and method that result in 
improved predictability of outcomes of corneal surgery. 
10 It is another object of the invention to provide an improved method and apparatus 

for design of optical surgery that minimizes invasiveness of the surgical procedure. 

It is another object of the invention to provide a method and apparatus for surgical 
design that results in reduction or elimination of postoperative irregular astigmatism. 

It is another object of the invention to provide an improved apparatus and method 
15 for surgical design which results in reduced multi-focal imaging of the central cornea, 
thereby enhancing contrast sensitivity and improving vision under low light illumination 
conditions. 

It is another object of the invention to provide an improved finite element analysis 
model of the human eye, including back-calculation of values of strain properties of the 
20 cornea and sclera, which incorporate the calculated strain properties of that eye and more 
accurately predict deformations of the cornea due to a hypothetical group of modeled 
incisions and/or excisions and/or ablation and/or insertions than has been achieved in the 
prior art. 
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It is another object of the invention to provide a system and method for providing 
an optimal surgical design for a human eye to achieve desired optical characteristics 
thereof. 

It is another object of the invention to reduce the likelihood of postoperative 
5 complications in the eye including, but not limited to over-correction or under-correction 
of pre-existing refractive errors. 

It is another object of the invention to provide a "training tool" or "surgery 
simulator" for surgeons who need to gain experience with corneal refractive surgery. 

It is another object of the invention to provide a device for designing new surgical 
10 procedures without the need for experimentation on live human beings. 

Briefly described, and in accordance with one embodiment thereof, the invention 
provides a system for simulating deformation of a cornea as a result of corneal incisions, 
excisions, ablations, and insertions in order to effectuate automated "surgical design" of a 
patient's eye in response to calculated strain conditions of the patient*s eye. A finite 
15 element analysis (FEA) model of the eye is constructed. Measured x,y,z coordinate data 
are interpolated and extrapolated to generate "nearest-fix" x,y,z, coordinates for the nodes 
of the finite element analysis mode. Measured thicknesses of the eye are assigned to each 
element of the finite element model. Pre-operative values of curvature of the cornea are 
computed. In one embodiment of the invention, strain property values are "back- 
20 computed" from measured stress values of corneal deformations at different pressure 
loads. An initial estimated surgical plan, including a number of incisions, locafions of 
incisions, incision orientations, incision depth, incision lengths, insert sizes, insert shapes, 
and insert locations is introduced into the shell finite element analysis model by 
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introducing duplicate "nodes" and nonlinear springs along the initial hypothetical 
incisions. Or, ablations may be included in the estimated surgical plan introduced into 
the finite element analysis model by varying the thickness and/or material property 
constants of the elements in the ablated region. A geometrically and materially nonlinear 
5 finite element analysis then is performed by solving the equations representing the finite 
element analysis model in response to incremental increases in intraocular pressure until 
the final "equilibrium state" is reached. Postoperative curvatures of the cornea are 
computed and compared to pre-operative values and to vision objectives. If the vision 
objectives are not met, the surgical model is modified and the analysis is repeated. This 

10 procedure is continued until the vision objectives are met. In one embodiment, a 
boundary element analysis model is used instead of a finite element analysis model. 

The novel features that are considered characteristic of the invention are set forth 
with particularity in the appended claims. The invention itself, how^ever, both as to its 
structure and its operation together with the additional object and advantages thereof v^ill 

15 best be understood from the following description of the preferred embodiment of the 
present invention when read in conjunction with the accompanying drawings. Unless 
specifically noted, it is intended that the words and phrases in the specification and 
claims be given the ordinary and accustomed meaning to those of ordinary skill in the 
applicable art or arts. If any other meaning is intended, the specification will specifically 

20 state that a special meaning is being applied to a word or phrase. Likewise, the use of the 
words "function" or "means" in the Description of Preferred Embodiments is not 
intended to indicate a desire to invoke the special provision of 35 U.S.C. §1 12, paragraph 
6 to define the invention. To the contrary, if the provisions of 35 U.S.C. §112, paragraph 
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6, are sought to be invoked to define the invention(s), the claims will specifically state the 
phrases "means for" or "step for" and a function, without also reciting in such phrases 
any structure, material, or act in support of the function. Even when the claims recite a 
"means for" or "step for" performing a function, if they also recite any structure, material 
5 or acts in support of that means of step, then the intention is not to invoke the provisions 
of 35 U.S.C. §112, paragraph 6. Moreover, even if the provisions of 35 U.S.C. §112, 
paragraph 6, are invoked to define the inventions, it is intended that the inventions not be 
limited only to the specific structure, material or acts that are described in the preferred 
embodiments, but in addition, include any and all structures, materials or acts that 
10 perform the claimed function, along with any and all knovra or later-developed 
equivalent structures, materials or acts for performing the claimed function. 

Brief Description of the Drawings 

FIG. 1 is a block diagram illustrating the components used in the invention. 
FIG. 2 is a basic flow chart useful in describing the method of the invention. 
15 FIG. 3 is a block diagram of a subroutine executed in the course of executing block 35 of 
FIG. 2 to interpolate and extrapolate data in order to obtain the nodal coordinates 
of a finite element analysis model. 
FIG. 4 is a block diagram of another subroutine executed in the course of executing block 
35 or FIG. 2 to "construct" the finite element analysis model. 
20 FIG. 5 is a three-dimensional diagram of the finite element mesh used in accordance with 
the present invention. 
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FIG. 6 is a partial side view illustrating both initial topography values of a portion of the 
cornea and final topography values resulting from simulated radial incisions and 
computed in accordance with the present invention. 

FIG. 7 is a diagram useful in explaining how incisions are included in the finite element 
5 analysis model of the present invention. 

FIG. 7A is a diagram useful in conjunction with FIG. 7 in explaining modeling of 
incisions. 

FIG. 8 is a diagram useful in explaining a technique for cubic spline interpolation and 
extrapolation to create "smoothed" three-dimensional data points from raw data 
1 0 provided by a keratoscope. 

FIG. 9 is a diagram useful in explaining automated back-calculation of the modulus of 
elasticity of the eye. 

FIG. 1 0 is a diagram useful in explaining optimization of the surgical plan according to 
block 41 of FIG. 2. 

15 Description of Preferred Embodiments 

The present invention involves constructing a strain determining model of a 
human eye using a suitable three-dimensional finite element analysis (FEA) model that 
includes a mesh that generally corresponds to the shape of the human eye. The finite 
element mesh is obtained using back calculated strain data and translated into the nodal 
20 points of the FEA model and describes the strain characteristics of the human eye. The 
nodal points in a small region are connected to each other, to form a finite set of 
elements. The elements are connected to each other by means of sharing common nodes. 
The strain values at any particular region are obtained by back calculation and are applied 
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to the elements. The "loading" of the finite element mesh structure is represented by the 
intraocular pressure, and the resistance of the structure to such applied "loading" is 
measured by the stiffness of the structure, which is computed on the basis of its 
geometry, boundary conditions, and its material properties, namely Poisson's ratio V, and 
5 Young's modulus E. 

In the area of structural mechanics, finite element analysis formulations are 
usually based on the "principle of virtual work," which is equivalent to invoking the 
stationary conditions of the total potential energy, D, given by 




equation 1 
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where 



s = BZ 



(2) 



and 



a=D8 



(3) 



8 is the transpose of the strain vector, Z is the transpose of nodal displacement 
15 vector, and Z^^ is the transpose of the nodal displacement vectors on the surface. Z and 
Z are nodal displacement terms associated with nodal loads. In the above equations, the 



various symbols have the following meanings: 



8 represents the strain vector 



D represents the material matrix 
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Z represents the vector of nodal displacements 



represents the nodal body force vector 



represents the nodal surface traction vector 
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dV represents differential volume 



dS represents differential surface area 



a represents the stress vector 



B represents a strain-displacement matrix 
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V represents volume 



S represents surface area. 

The first term on the right hand side of the equation (1) is the strain energy of the 
structure, and the second and third terms represent the total work accomplished by the 
external forces and body forces. The strain energy is a function of the strains and stresses 
10 that are related to each other via the material matrix D. The material properties that 
contribute to the material matrix D include the modulus of elasticity (Young's modulus) 
and Poisson*s ratio. In a uniaxial state of stress, Poisson's ratio is defined as: 

eiat =-VSiong, (4) 

where siat is the normal strain in the lateral direction and siong is the normal strain 
15 in the longitudinal direction. In a uniaxial state of stress, Young*s modulus E is defined 
according to 

CTZZ =E8zz, (5) 

where azz is the normal stress and Szz is the normal strain. The work 
accomplished is a function of the applied loads and surface tractions. 
20 Using an assumed displacement field, the minimization of the total potential 

energy n leads to the element equilibrium equations of the form 




(6) 



where the expression of 
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nxn 



(l/v}B'^DBdV 



is the element stiffness matrix, and Znxi is the vector of element nodal 
displacements rnxi is the vector of element nodal forces. Since the entire structure is 
assumed to be in equilibrium, the assembly of the element equations leads to the 
5 structural equilibrium equations of the form 



where Knxn is the structural stiffness matrix, Znxi is the vector of nodal 
displacements and Rnxi is the vector of nodal forces. These algebraic equations are 
finally solved for Z in a variety of ways depending on whether the structural behavior is 

1 0 linear or nonlinear. 

A commercially available finite element analysis program that effectively solves 
these equations after the appropriate values and boundary conditions have been assigned 
to the various nodes and the appropriate material properties have been assigned to the 
various elements defined by the connectivity of the nodes is called ABAQUS, available 

15 from HKS, Inc. of Providence, R.I. Creating the FEA model for purposes of the present 
invention simply involves inputting to the ABAQUS program the x, y, z coordinates for 
each node, inputting the strains that act on the nodes and/or elements, assigning 
appropriate boundary conditions to each node, defining the nodal connectivity that 
defines each element, and inputting the eye material properties and thickness or stiffness 

20 to each defined element along with other input data, such as whether the analysis is linear 
or non-linear, or the properties and definitions of the non-linear springs. 

It should be noted that there are two popular approaches to solving finite element 
analysis problems, one being the above-described approach of minimizing total potential 



KnxNXZnxI -RnxI 



(7) 
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energy (or, the variational approach), the other being a method of weighted residuals 
which operates on partial differential equations defining the problem. The first approach 
is generally recognized to be simpler, and is implemented by the above ABAQUS 
program, but the invention could be implemented using the second approach. 
5 FIG. 1 shows an apparatus used in conjunction with the present invention. An 

ultrasonic instrument 15, such as a DGH packymetor, model DGH-2000 available from 
DGH Technology, Inc. of Frazier, Pa., is used to obtain the thickness and intraocular 
pressure of cornea 11 A. 

A corneal topographer 12, which can be a model TMS-1, manufactured by 

10 Computed Anatomy, 28 West 36th Street, New York, N.Y., is utilized to measure the 
surface topography of cornea 11 A. The resulting information is transferred by means of 
a floppy disk to a computer system 14. Alternatively, a digital data bus 13 could be 
provided to transfer topography information from corneal topographer 12 to computer 
system 14. A printer 17 is connected by a cable to the printer port of the computer 

15 system 14. 

Thickness and interocular pressure measurements are made by ultrasonic 
instruments 15. This data then is used in the generation of the finite element model. 
However, it is possible to have this data transferred digitally, either by means of a floppy 
disk or a communication link, to the personal computer 14. 
20 A conventional pressure loading device 19 is utilized to apply a precisely 

measurable force on a point of the sclera as far away as practical from the cornea, so that 
resulting changes on the elastic cornea as a result of the new loading can be measured. 
Then, in accordance with the present invention, the value of Young's modulus can be 
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"back-calculated" in the manner subsequently described. Alternately, uniform pressure 
loading could be achieved by applying a sealed pressure chamber to the eye and 
increasing the gas pressure therein. Such uniform loading may have the advantage of 
providing less "noise" error in the measurements. A suitable pressure loading device 19 
5 could be an opthalmo dynamometer, commercially available from Bailliart, of Germany. 

To obtain an FEA model of the patient*s eye, the measured topographical data is 
interpolated and extrapolated using the subsequently described cubic spline technique to 
provide a pre-established reduced number of nodal points of a finite element mesh, with 
nodal coordinates which are a "close fit" to the measured corneal surface. Values of the 

10 thickness of the cornea and sclera obtained from the data obtained from ultrasonic 
instrument 15 are assigned to the various finite elements of the FEA model. The FEA 
mesh then defines a continuous surface that accurately represents the pre-operative 
surface of the cornea, including any astigmatism that may be present. 

The curvatures of the surface then are computed at each node of the finite element 

15 analysis model. Surfaces of revolution are generated by revolving a plane curve, called 
the meridian, about an axis not necessarily intersecting the meridian. The meridian 
(defined by a radial line such as 21 in FIG. 5) is one of the principal sections and its 
curvature at any point is one of the principal curvatures kj. (The principal curvatures are 
defined as the maximum and the minimum curvatures at a point on a surface.) The other 

20 orthogonal principal section is obtained by the intersection of the surface with a plane 
that is at right angles to the plane of the meridian and that also contains the normal. The 
second principal section has curvature k2. If the equation of the meridian is written as r = 
f(z), then 
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K, =-r"/[l+(r')Y^ 
and 

\^2 = [r[\Hf)fY 

r' and r" being the first and second derivatives of r, respectively. 
5 FIG. 2 is a flowchart useful in explaining the basic steps involved in use of the 

system shown in FIG. 1 to produce an optimum design for guiding surgery of a patient's 
eye. In block 31 of FIG. 2, the physician determines the "vision objectives" for the eye. 
The vision objectives can be specified as homogenous strain relationships throughout the 
cornea when the eye is in an accommodatively relaxed condition after completion of the 

10 surgery. The strains at various locations can be displayed, for example, by graded colors. 
These desired strain changes may be determined at nodal points on the cornea by back 
calculation and use of a spatially resolved refi'actometer. The vision objectives are 
selected to maximize the number of light rays that the eye focuses on the fovea for a 
given functional distance by creating a homogenous or uniform strain field in the cornea. 

15 In block 32, the physician provides initial estimates of the number of incisions, 

ablations, insertions or thermal shrinkages required, their locations, the orientations of the 
various incisions, ablations, insertions or thermal shrinkages, the incision, ablation, or 
thermal shrinkage lengths, the incision, ablation, or thermal shrinkage depths, and the 
size and shape of the insertions needed in order to accomplish the vision objectives of 

20 block 31. As indicated in block 33, a corneal topographer 12 is used to obtain a 
topographic map of a portion of the eye. The TMS-1 corneal topographer mentioned 
above is capable of providing an x, y, z coordinate "map" that covers most of the cornea. 
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producing a data file from which the x, y, z coordinates of approximately 7000 points can 
be obtained. 

As indicated in block 34 of FIG. 2, the ultrasonic instrument 15 is used to provide 
measurements of the thickness of the cornea and the intraocular pressure. In the 
5 prototype system presently being implemented, typical values of Poisson's ratio and 
Young's modulus are used. Presently, Poisson's ratio values of 0.49 are used for both the 
cornea and sclera. Presently, values of Young's modulus equal to 2 10^ dynes per square 
millimeter are used for the cornea and 5 10^ dynes per square millimeter for the sclera. 

Preferably, Young's modulus, and ultimately elemental strain, is "back- 

10 calculated" on the basis of corneal topographical changes measured by using the corneal 
topographer 12 after varying a known force applied by pressure loading device 19 (FIG. 
1) to the eye. The main objective of the back-calculation procedure is to determine as 
accurately as possible the modulus of elasticity for the cornea and the sclera, because it 
also is recognized that these values vary from patient to patient, and because it also is 

15 recognized that the modulus of elasticity is one of the most crucial parameters that 
influences the finite element analysis predictions. To describe the basic technique of the 
back-calculation procedure, refer to FIG. 9, which shows three assumed states, namely 
State 0 in which the cornea is relaxed. State 1 in which pressure loading device 19 applies 
point load PI to the sclera of eye 11, and State 2 in which pressure loading device 19 

20 applies point load P2 to the sclera. PO is the intraocular pressure that is uniformly applied 
to the inner surface of the cornea and sclera. The values of the moduli of elasticity for 
the cornea and the sclera, respectively, are adjusted such that the z coordinates at selected 
nodes are close to the values actually measured for State 1 and State 2 by TMS-1 corneal 
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topographer 12 with the two values of point loads PI and P2 actually applied, 
respectively. 

Let Zij be the observed z coordinate at a particular node i for a state j obtained 
using the above mentioned TMS-1 system. Let Zy be the computed z coordinate at a 
5 particular node i for the state j using the finite element analysis according to the present 
invention. The back-calculation problem then is to find the value of {Ec, Es} to minimize 
the value of the expression for 



Es<'^Es<Es^ 

and where {Ec,Es} is the vector of design variables, f(Ec,Es) is the objective 
function, Ec is the modulus of elasticity of the cornea, Es is the modulus of elasticity of 

1 5 the sclera, n is the number of points at which z displacements are to be computed, and the 
two inequality constraints represent the lower (L) and upper (U) bounds on the two 
design parameters. It should be appreciated that such a problem formulation falls under 
the category of a non-linear programming (NLP) problem, and can be solved using 
various non-linear programming techniques such as the "method of feasible directions", 

20 or using a constrained least-squares technique. A commercially available program for 
solving such non-linear programming problems is the DOT (Design Optimization Tools) 
program, available from VMA, Inc. of Santa Barbara, Calif. 



f(Ec,Es) = 5:S(-l +Zf'/ZJ)2 



with the conditions 



10 



Ec^<Ec<Ec" 



u 



and 
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As indicated in block 35, an FEA model is "constructed" for the patient's eye by 
interpolating between the various 7000 x, y, z coordinates of the corneal map produced 
by TMS-1 corneal topographer 12 to provide a smaller number of representative 
"smoothed" x, y, z values to be assigned to the various nodes of the FEA mesh shown in 



FIG. 5 shows one quadrant of the FEA mesh, the other three quadrants being 
substantially identical except for the nodal values assigned to the nodes thereof. The 
FEA mesh shown in FIG. 5 includes a plurality of equi -angularly spaced radial lines 21, 
each extending from a cornea center or apex 27 of cornea section 24 to the bottom of 

10 sclera section 23. In the FEA mesh actually used in a prototype of the invention under 
development, there are 32 such radial lines 21 and also 30 generally equally spaced 
circumferential lines 22. Each area such as 25 that is bounded by two adjacent radial 
lines and two adjacent circumferential lines 22 is an "element" of the FEA model. Each 
typical element 25 includes eight assigned "nodes", such as nodes 26-1 ... 26-8. The 

1 5 four comers of a typical element such as 25 share comer nodes 26-1 , 3, 5, 7 with adjacent 
elements, and also share "midpoint" nodes 26-2, 4, 6, and 8 with corresponding midpoint 
nodes of adjacent elements. The nodes and the cormectivity thereof which define the 
elements of the FEA mesh thus are illustrated in FIG. 5. 



20 coordinates and its boundary conditions, which are whether the node can or cannot 
undergo x, y, z displacements and rotations. The values assigned to each element in the 
FEA model include the thickness of the element, Young's modulus or modulus of 
elasticity, the shear modulus, the strain, and Poisson*s ratio in the orthotropic directions, 



5 



FIG. 5. 



The values assigned to each node include its interpolated/extrapolated x, y, z 
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namely the xy, xz, and zx directions. Any external "loading" forces at each node also are 
assigned to that node. The orthotropic values of Poisson's ratio presently uses are Uxy 
=0.0025, Uxz =0.0025, and Uzx =0.49. The value of shear modulus used is Gyz =6.71 10^ 
dynes per square millimeter. 
5 The objective of the tasks in block 35 is three-fold. First, the total number of 

nodes, and thus the elements generated from them, should be a variable, so that the mesh 
sensitivity of the results can be studied while the operator is given the chance to use a 
coarse mesh for preliminary studies. Second, the nodal points generated should be 
compatible With the choice of element required. For example, eight-node shell elements 

10 are used in the present approach. However, the proposed system is able to generate any 
type of element required, such as a 27-node hexahedral three-dimensional element, a 6- 
node triangular shell element, or a 9-node shell element. Third, the nodes generated must 
be able to provide sufficient mesh refinement or density to achieve the needed resolution. 
A refined mesh in the regions of primary interest such as the optical zone is important 

15 since it can capture the stresses, strains and the variations of the displacements, and thus, 
the curvatures. The mesh refinement parameter is chosen by the operator as one of the 
variables to study for the regions of primary interest, such as the optical zone (i.e., the 
portion of the cornea central to the radial incisions), while the other regions such as the 
sclera are still incorporated in the model. 

20 Some of the steps performed by computer 14 in accordance with block 35 of FIG. 

2 are shown more specifically in FIGS. 3 and 4. As indicated in block 45 of FIG. 3, 
computer 14 reads the ASCII data files containing the above-mentioned 7000 coordinates 
of the corneal map produced by the TMS-1 comeal topographer 12. Due to the nature of 
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the data collection, it is possible that some "noise" exists in the original data. The origin 
of the noise might be attributed to the inability of corneal topographer 12 to provide an 
exact determination of the coordinates, or the lack of existence of the coordinate value at 
an expected site. As indicated in block 46, a simple program scans the original data for 
5 elimination of such data points. 

As indicated in block 47 of FIG. 3, the polar coordinate data supplied by the 
TMS-1 corneal topographer is converted into the above-mentioned 7000 x, y, z 
coordinates. Points which lie along the radial lines 21 of the FEA mesh shovvn in FIG. 5 
are selected for use in the interpolation/extrapolation process described below. 

10 As indicated in block 48, the cubic spline interpolation and extrapolation 

procedure (described later with reference to the diagram of FIG. 8) is utilized to compute 
the intermediate x, y, z coordinates for each node of the FEA mesh lying on the pre- 
defined radial line 21 (FIG. 5). Then, as indicated in block 49 of FIG. 3, the program 
creates a final set of x, y, z coordinates for nodes that lie on the circumferential lines 22 

1 5 of the FEA mesh (FIG. 5) using the cubic spline interpolation/extrapolation method. This 
step is necessary since data points that have the same radial coordinate do not necessarily 
have the same height or z value. 

At this stage, as indicated in block 50 of FIG. 3, the initial or pre-surgery diopter 
values at the final setup points are computed. Once the radial lines 21 are "generated", a 

20 series of nodes are selected at a specific height and used to obtain the circumferential 
nodes, i.e., the nodes which are on circumferential lines 22 "between" the radial lines 21. 
The x, y, z coordinates and diopter values of curvature at each node of the model then are 
output. 



21 



As indicated in block 52 the coordinate data files produced by block 51 of FIG. 3 
are read. In block 53 the finite element mesh options/data are read. This includes the 
orthotropic material properties of the cornea and sclera, the nonlinear load-elongation 
curve data used by the "spring" elements (i.e., the insertion thickness or thermal 
5 shrinkage depth as subsequently described with reference to FIGS. 7 and 7 A), the loading 
information (i.e., the intraocular pressure), and the boundary conditions (i.e., the 
connections of the bottom nodes of the sclera to a stationary reference). Then, the 
program reads the surgical data, as indicated in block 54, and goes to block 55 in which 
the FEA model is "created", i.e., duplicate nodes, element connectivity and load- 

1 0 elongation data for the spring elements are created. 

Finally, in block 56, the data files required for carrying out a geometric and 
materially nonlinear finite element analysis are created and output. In the present 
embodiment of the invention, the above-mentioned ABAQUS program is used as the 
finite element analysis program and is executed on the computer system 14. 

15 Returning to FIG. 2, in block 36, pre-operative curvatures are completed in 

diopters at each node of the FEA model. 

Then, in block 37, an initial (estimated) number of insertions or thermal 
shrinkages are "constructed" in the FEA model using the information established in block 
32. It should be recognized that incisions and ablations or linear combinations of the 

20 entire above are included herein, but for simplicity, the following discussion will be 
limited to insertions and thermal shrinkages. FIGS. 7 and 7A illustrate how each such 
thermal shrinkage is modeled in accordance with the present invention. In FIG. 7, the 
FEA mesh 11 A includes a thermal shrinkage 61 modeled along a radial line 58 of the 
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FEA mesh, in which numerals 61-1,2,3,4,5 represent all nodes of the FEA mesh from one 
end of the modeled thermal shrinkage 61 to the other. The elasticity at these nodes, and 
their neighboring spring elements, are changed according to predetermined values 
corresponding to the changes caused by thermally induced shrinkage of the collagen 
5 fibers of the cornea. When thermal shrinkage of the collagen causes a stiffening of the 
comeal tissue, the elasticity is reduce; commensurately, when thermal shrinkage of the 
collagen causes a "loosening" of the comeal tissue, the elasticity is increased. The 
thermally induced shrinkages may be caused by laser heating, cauterizing wires, or the 
like. 

10 These changed spring elements have nonlinear load-deflection curves. The nature 

of the curves is a function of the depth of thermal shrinkages and the material properties 
of the tissue through which the thermal shrinkage is made. The depth of the modeled 
thermal shrinkage 61 is represented by equations corresponding to the nonlinear elastic 
spring elements in FIG. 7A. When an FEA program is executed, the effect of the 

15 intraocular pressure is to cause the thermal shrinkage 61 to change an amount determined 
by the elastic spring constants assigned to the neighboring spring elements. 

In another example FIGS. 7B and 7C illustrate how each such insertion is 
modeled in accordance with the present invention. In FIG. 7B, the FEA mesh 11 A 
includes an insertion 61B modeled along several radial lines of the FEA mesh, in which 

20 numerals represent all nodes of the FEA mesh impacted by the 

modeled insertion 61B. In this example, z value for the impacted nodes and the elastic 
constants are at the nodes, including their neighboring spring elements, are changed by 
the addition of the insert into the comeal tissue. The inclusion of the insert effectively 
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stiffens the neighboring spring elements, similar to the effect of increasing pressure. The 
size, depth, and shape of the insert may be varied, as desired, and may be preferably 
designed to ultimately provide for a homogeneous strain relationship in the corneal 
tissue. 

5 These changed spring elements have nonlinear load-deflection curves. The nature 

of the curves is a function of the size, depth, and shape of the insert and the material 
properties of the tissue through which the thermal shrinkage is made. The size, depth, 
and shape of the modeled insert 61B is represented by equations corresponding to the 
nonlinear elastic spring elements in FIG. 7C. When an FEA program is executed, the 
10 effect of the intraocular pressure is to that caused by the insert 61B. 

The above-mentioned ABAQUS FEA program, when executed as indicated in 
block 37 of FIG. 2, computes the displacements at each node of the FEA model in 
response to the intraocular pressure. 

The computed nodal x, y, z displacements are added to the corresponding pre- 
15 operative x, y, z values for each node, and the results are stored in a data file. If desired, 
the results can be displayed in, for example, the form illustrated in FIG. 6. Post-operative 
curvatures (computed in diopters) and corneal strains then are computed and displayed 
for each node based on the new nodal locations. 

In FIG. 6, which shows a computer printout produced by the system of FIG. 1, the 
20 measured pre-operative configuration of the eye surface is indicated by radial lines 82, 
and the computed post-operative configuration is indicated by radial lines 84. Numerals 
21 generally indicate radial lines of the FEA model, as in FIG. 5. More specifically, 
numerals 21A-1 and 21 A-2 designate radial lines of the pre-operative surface represented 
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in the FEA model, and numerals 21B-1 and 21B-2 represent radial lines of the 
"computed" post-operative surface in the FEA model. Numerals 61 A designate proposed 
radial incisions in the "measured" pre-operative surface, and numerals 61 B designate the 
same incisions in the "computed" post-operative surface. (The individual circumferential 
5 lines of the computer printout of FIG. 6 are difficult to identify, but this does not prevent 
accurate interpretation of the effect of the proposed incisions on the curvature of the 
cornea.) Numeral 86 indicates the limbus. 

As indicated in block 38, the post-operative curvatures and strains are then 
compared with the pre-operative curvatures and strain with their corresponding vision 

10 objectives established according to block 31 to determine whether the initial estimated 
surgical plan accomplished the vision objectives. 

Then, as indicated in block 39, computer 14 determines if the strain boundary 
conditions along with the vision objectives are met. If the determination of decision 
block 39 is affirmative, the surgical design is complete, as indicated in label 40. 

15 Otherwise, however, the program executed by computer 14 goes to block 41, and an 
optimization technique, subsequently described with reference to FIG. 10, is utilized to 
modify the number of incisions, ablatoins, thermal shrinkages, and insert, their locations, 
orientations, lengths, depths, sizes, and shapes. The process then returns to block 37 and 
repeats until an affirmative determination is reached in decision block 39. 

20 The technique for modifying and optimizing the surgical design according to 

block 41 can be understood with reference to FIG. 10. As indicated above, the vision 
objectives are to obtain prescribed curvature values at specific FEA nodal locations i on 
the cornea. As an example, assume the surgical plan includes the locations of the 
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insertions and includes the location, size and shape of each insertion. The surgical 
optimization problem of block 41 then can be defined to be the problem of determining a 
value of {aj,lj,dj} that minimizes the value of the expression 

F(aj,lj,d,) = Z(-l+ri/ri) 
5 with 

aj^ < aj < Aj^ 

Ij^ < Lj < Lj^ 

dj^<dj<dj^ 

where {aj,lj,dj} is the vector of design variables, f{aj,lj,dj} is the objective 
10 function, aj is the starting radial distance from the center of the finite element model, as 
shown in FIG. 10, Ij is the length of the insertion, dj is the height of the insertion, j is the 
insertion number, ri is the computed curvature value based on the results fi-om the finite 
element analysis, ri is the observed curvature, n is the number of points at which 
curvature computations are to be carried out, and the three above inequality constraints 
1 5 represent the lower (L) and upper (U) bounds on the three design parameters. 

With this information, it is possible to include any parameter that influences the 
finite element model as a potential design variable, and any response or parameters 
related to the response computed by the finite element analysis as appearing in the 
objective functions or constraints. For example, the shape of the insertion, the number of 
20 different insertions, or the compressibility of a thickness of tissue can be design variables. 

The foregoing problem formulation falls under the category of nonlinear 
programming problem. Those skilled in the art can readily solve such problems using 
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nonlinear programming techniques utilizing commercially available nonlinear 
programming software, such as the previously mentioned DOT program. 

FIG. 8 is usefiil in illustrating the application of cubic spline techniques referred 
to in block 48 of FIG. 2 to interpolate/extrapolate data from the nodal points of the FEA 
5 mesh from the data points obtained from TMS-1 corneal topographer 12. In FIG. 8, 
numeral 71 designates the z axis or a center line of the cornea passing through its apex, 
numeral 72 designates a radial line along which data points obtained from corneal 
topographer 12 lie, and numeral 73 designates various such measured data points. The 
extent of the cornea is indicated by arrow 78, and the extent of the sclera is indicated by 
10 arrow 79. The extent of the "optical zone" is indicated by arrow 80. As indicated above, 
the TMS-1 corneal topographer provides 7000 such data points 73. The first step of the 
cubic spline process takes such data points, as indicated by arrow 74, and "fits" each 
segment of radial line 72 between adjacent corneal topographer data points 73 to the 
equation 

15 z=ax^ +bx^ +cx+d, 12 ^ 

where z is a distance along center line 71, and x is distance in the horizontal 
direction from line 71 toward the base of the sclera. Equation 12 then is used to compute 
values of z for each value of x corresponding to a node of the FEA mesh (shown in FIG. 
5) to obtain values z for each of the nodes of the FEA mesh along each radial line 21, as 

20 indicated by arrow 76 in FIG. 8. Values for the "midpoint" nodes such as 26-2 and 26-6 
of FIG. 5 are obtained by interpolating adjacent nodal values of z on the same 
circumferential line 22. Most texts on numerical analysis disclose details on how to use 
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the cubic spline technique, and various commercially available programs, such as IMSL, 
available from IMSL, Inc. of Houston, Tex., can be used. 

A fixed boundary condition for the base of the sclera can be assigned. It has been 
found that the nature of the boundary conditions at the base of the sclera has only a small 
5 effect on the results of the finite element analysis of the cornea. 

The above-described model was utilized to compute the strains and nodal 
deflections in a particular patient's eye based on measured topographical data extending 
outward approximately 8 millimeters from the center of a patient*s eye. The measured 
data was extrapolated outward another 8 millimeters to approximate the topography of 

1 0 the remaining cornea. 

The above-described FEA model can be used to pre-operatively design incisions, 
excisions or ablations, thermal shrinkages, and insertions into the cornea, resulting in 
great predictability of surgical outcome and thereby allowing minimum invasiveness to 
achieve the desired result with the least amount of surgical trauma to the cornea. Fewer 

15 operative and post-operative visits by the patient to the surgery clinic are likely as a result 
of the use of this procedure. Advantages of the improved surgical designs that result 
from the above-described invention include reduced multi-focal imaging of the central 
cornea, thereby enhancing contrast sensitivity and improving vision under low light 
illumination conditions. Reduction or elimination of post-operative irregular astigmatism 

20 is another benefit. Yet another benefit is minimization of side effects such as glare and 
fluctuation of vision associated with traditional incisional keratotomy. The described 
mathematical model will have other uses, such as allowing design of a bifocal corneal 
curvature to allow both near and distance vision for patients in the presbyopic stage of 
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their lives. The model of the present invention also will allow development of new 
surgical techniques for correcting nearsightedness, farsightedness and astigmatism as a 
viable alternative to experimenting on live human corneas. 

While the invention has been described with reference to several particular 
5 embodiments thereof, those skilled in the art will be able to make the various 
modifications to the described embodiments of the invention without departing from the 
true spirit and scope of the invention. It is intended that all combinations of elements and 
steps which perform substantially the same function in substantially the same way to 
achieve the same result are within the scope of the invention. 

10 For example, keratoscopes or other cornea measurement devices than the TMS-1 

device can be used. Non-radial incisions, such as T-shaped incisions for correcting 
astigmatism, can be readily modeled. Many variations of the finite element model are 
possible. In the two-dimensional shell finite element analysis model described above, the 
use of the nonlinear springs to model depths of incisions could be avoided by modeling 

15 elements around the proposed incision to have reduced thickness and/or different material 
properties, so that the incision region has reduced stiffness, and the computed 
deformations are essentially the same as if the nonlinear springs were to be used. For 
example, it is possible to use three-dimensional finite elements in lieu of the two- 
dimensional shell finite elements with assigned thickness parameters, and model the 

20 incisions directly, without having to use the nonlinear spring elements. Mathematical 
models other than a finite element analysis model can be used. For example, a boundary 
element analysis model could be used. As those skilled in the art know, the basic steps in 
the boundary element methods are very similar to those in the finite element methods. 
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However, there are some basic differences. First, only the boundary is discretized, that 
is, the elements are "created" only on the boundary of the model, whereas in finite 
element analysis models the elements are "created" throughout the domain of the model. 
Second, the fundamental solution is used which satisfies the governing differential 
5 equation exactly. A fundamental solution is a function that satisfies the differential 
equation with zero right hand side (i.e., with body force set to zero) at every point of an 
infinite domain except at one point knovm as the source or load point at which the right 
hand side of the equation is infinite. Third, the solution in the interior of the model can 
be obtained selectively once the approximate solution on the boundary is computed. 

10 Although constant intraocular pressure has been assumed, non-constant intraocular 
pressure could be incorporated into the described technique. Although post-operative 
swelling has been assumed to not effect the eventual curvatures of the cornea, healing of 
the incision does effect the eventual curvature. The finite element analysis model can be 
adapted to model such healing effects and predict the final curvatures, strains, etc. 

15 Additionally, p-finite elements, Raleigh-Ritz, mixed formulations, Reissner*s 

Principal, all can be used to generate the finite element equations. These equations, then, 
can be used in the modeling method of the present invention. 

The preferred embodiment of the invention is described above in the Drawings 
and Description of Preferred Embodiments. While these descriptions directly describe 

20 the above embodiments, it is understood that those skilled in the art may conceive 
modifications and/or variations to the specific embodiments shown and described herein. 
Any such modifications or variations that fall within the purview of this description are 
intended to be included therein as well. Unless specifically noted, it is the intention of 
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the inventor that the words and phrases in the specification and claims be given the 
ordinary and accustomed meanings to those of ordinary skill in the applicable art(s). The 
foregoing description of a preferred embodiment and best mode of the invention known 
to the applicant at the time of filing the application has been presented and is intended for 

5 the purposes of illustration and description. It is not intended to be exhaustive or to limit 
the invention to the precise form disclosed, and many modifications and variations are 
possible in the light of the above teachings. The embodiment was chosen and described 
in order to best explain the principles of the invention and its practical application and to 
enable others skilled in the art to best utilize the invention in various embodiments and 

10 with various modifications as are suited to the particular use contemplated. 
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